{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ch 5.6 改进的BP算法\n",
    "试设计一个BP改进算法，能通过动态调整学习率显著提升收敛速度。编程实现该算法，并选择两个UCI数据集与标准BP算法进行比较。\n",
    "\n",
    "Note：\n",
    "1. 学习率调整方法包括Adagrad，RMSProp等\n",
    "2. 这里使用AdaGrad算法，算法如下（Ian Goodfellow《深度学习》P188）：\n",
    "常数sigma=1e-7\n",
    "\n",
    "梯度累积变量 r = 0\n",
    "\n",
    "全局学习率 epsilon = 0.1\n",
    "\n",
    "while 没有达到停止准则 do\n",
    "\n",
    "    采集m个样本\n",
    "    \n",
    "    计算m个样本的平均梯度grad\n",
    "    \n",
    "    累积平方梯度r = r + grad * grad\n",
    "    \n",
    "    计算更新theta = theta - epsilon/(sigma+sqrt(r)) * grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dataset = pd.read_csv('data/blood.txt',header=None)\n",
    "dataset.columns = ['x1','x2','x3','x4','label']\n",
    "dataset = dataset.apply(lambda x:(x-np.min(x))/(np.max(x)-np.min(x)))\n",
    "dataset_array = np.array(dataset)\n",
    "dataset_array = dataset_array[np.random.permutation(len(dataset))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    return 1/(1+np.exp(-x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def cal_acc(dataset_array, w_h_j, v_i_h, theta_j, gamma_h):\n",
    "    dataset = dataset_array[:,:-1]\n",
    "    label = dataset_array[:,-1]\n",
    "    acc = 0\n",
    "    for i in range(len(dataset)):\n",
    "        pred_label = predict(dataset[i], w_h_j, v_i_h, theta_j, gamma_h)\n",
    "        if pred_label == label[i]:\n",
    "            acc += 1\n",
    "    return acc/len(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def predict(x_i, w_h_j, v_i_h, theta_j, gamma_h):\n",
    "    alpha_h = np.dot(x_i, v_i_h)\n",
    "    b_h = sigmoid(alpha_h - gamma_h)\n",
    "    beta_j = np.dot(b_h, w_h_j)        \n",
    "    y_j_cap = sigmoid(beta_j - theta_j)\n",
    "#     print(y_j_cap, end=\" \")\n",
    "    if y_j_cap[0][0] > y_j_cap[0][1]:\n",
    "#         print(0)\n",
    "        return 0\n",
    "    else:\n",
    "#         print(1)\n",
    "        return 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def adagrad_bp(dataset_array, niter):       \n",
    "    h_num = 30\n",
    "    j_num = 2\n",
    "    i_num = dataset_array.shape[1] - 1\n",
    "    l_num = 2\n",
    "    epsilon = 0.1\n",
    "    sigma = 10e-7\n",
    "    m_num = 20\n",
    "    w_h_j = np.random.random([h_num, j_num])\n",
    "    v_i_h = np.random.random([i_num, h_num])\n",
    "    theta_j = np.random.random([1,j_num])\n",
    "    gamma_h = np.random.random([1,h_num])\n",
    "    r_w_h_j = np.zeros([h_num, j_num])\n",
    "    r_v_i_h = np.zeros([i_num, h_num])\n",
    "    r_theta_j = np.zeros([1,j_num])\n",
    "    r_gamma_h = np.zeros([1,h_num])\n",
    "    error_list = []\n",
    "    for i in range(niter):\n",
    "        error = 0\n",
    "        delta_w_h_j = np.zeros([h_num, j_num])\n",
    "        delta_theta_j = np.zeros([1, j_num])\n",
    "        delta_v_i_h = np.zeros([i_num, h_num])\n",
    "        delta_gamma_h = np.zeros([1, h_num])\n",
    "        for _ in range(m_num):\n",
    "            k = np.random.randint(len(dataset_array))\n",
    "            x_i, y_i = dataset_array[k,:-1], dataset_array[k,-1]  # x_i：1xi的向量\n",
    "            x_i = np.array(x_i)\n",
    "            x_i = np.reshape(x_i, [1,i_num])\n",
    "            y_j = np.zeros([1,l_num])\n",
    "            y_j[0][int(y_i)] = 1\n",
    "            alpha_h = np.dot(x_i, v_i_h)\n",
    "            b_h = sigmoid(alpha_h - gamma_h)\n",
    "            beta_j = np.dot(b_h, w_h_j)        \n",
    "            # formula 5.3\n",
    "            y_j_cap = sigmoid(beta_j - theta_j)\n",
    "            # formula 5.10\n",
    "            g_j = y_j_cap * (np.ones_like(y_j_cap)-y_j_cap)*(y_j - y_j_cap)\n",
    "            # formula 5.15\n",
    "            e_h = b_h *(np.ones_like(b_h)-b_h) * np.dot(g_j, w_h_j.T)\n",
    "            \n",
    "            delta_w_h_j += np.dot(b_h.T, g_j)\n",
    "            delta_theta_j += -g_j\n",
    "            delta_v_i_h += (np.dot(e_h.T, x_i)).T\n",
    "            delta_gamma_h += -e_h\n",
    "            error += 0.5 * np.sum((y_j_cap - y_j) * (y_j_cap - y_j))\n",
    "\n",
    "        delta_w_h_j /= m_num\n",
    "        delta_theta_j /= m_num\n",
    "        delta_v_i_h /= m_num\n",
    "        delta_gamma_h /= m_num\n",
    "\n",
    "        r_w_h_j += delta_w_h_j*delta_w_h_j\n",
    "        r_v_i_h += delta_v_i_h*delta_v_i_h\n",
    "        r_theta_j += delta_theta_j*delta_theta_j\n",
    "        r_gamma_h += delta_gamma_h*delta_gamma_h\n",
    "\n",
    "        w_alpha = np.ones_like(delta_w_h_j) * epsilon / (np.ones_like(delta_w_h_j) * sigma + np.sqrt(r_w_h_j))\n",
    "        theta_alpha = np.ones_like(delta_theta_j)*epsilon/(np.ones_like(delta_theta_j)*sigma+np.sqrt(r_theta_j))\n",
    "        v_alpha = np.ones_like(delta_v_i_h)*epsilon/(np.ones_like(delta_v_i_h)*sigma+np.sqrt(r_v_i_h))\n",
    "        gamma_alpha = np.ones_like(delta_gamma_h)*epsilon/(np.ones_like(delta_gamma_h)*sigma+np.sqrt(r_gamma_h))\n",
    "\n",
    "        #这里用+是因为delta里面其实是负梯度\n",
    "        w_h_j += w_alpha*delta_w_h_j\n",
    "        theta_j += theta_alpha*delta_theta_j\n",
    "        v_i_h += v_alpha*delta_v_i_h\n",
    "        gamma_h += gamma_alpha*delta_gamma_h\n",
    "\n",
    "\n",
    "        if i % 10 == 0:\n",
    "            error_list.append(error/m_num)\n",
    "    print(error_list)\n",
    "    return w_h_j, v_i_h, theta_j, gamma_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.49580279793352455, 0.30749959061012194, 0.23420224929534378, 0.15867765393551489, 0.18506877978027941, 0.18003805913516061, 0.13009109482011777, 0.15181505713283036, 0.34733435432216569, 0.12839970876641127]\n",
      "0.8120805369127517\n",
      "Wall time: 135 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "test_num = dataset.shape[0]//5\n",
    "w_h_j, v_i_h, theta_j, gamma_h = adagrad_bp(dataset_array[test_num:], niter =100)\n",
    "acc = cal_acc(dataset_array[:test_num], w_h_j, v_i_h, theta_j, gamma_h)\n",
    "print(acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def standard_bp(dataset_array,niter =350, eta=0.1):    \n",
    "    h_num = 10\n",
    "    j_num = 2\n",
    "    i_num = dataset_array.shape[1] - 1\n",
    "    l_num = 2\n",
    "    w_h_j = np.random.random([h_num, j_num])\n",
    "    v_i_h = np.random.random([i_num, h_num])\n",
    "    theta_j = np.random.random([1,j_num])\n",
    "    gamma_h = np.random.random([1,h_num])\n",
    "    error_list = []\n",
    "    for i in range(niter):\n",
    "        error = 0\n",
    "        for k in range(len(dataset_array)):\n",
    "            x_i, y_i = dataset_array[k,:-1], dataset_array[k,-1]  # x_i：1xi的向量\n",
    "            x_i = np.array(x_i)\n",
    "            x_i = np.reshape(x_i, [1,i_num])\n",
    "            y_j = np.zeros([1,l_num])\n",
    "            y_j[0][int(y_i)] = 1\n",
    "            alpha_h = np.dot(x_i, v_i_h)\n",
    "            b_h = sigmoid(alpha_h - gamma_h)\n",
    "            beta_j = np.dot(b_h, w_h_j)        \n",
    "            # formula 5.3\n",
    "            y_j_cap = sigmoid(beta_j - theta_j)\n",
    "            # formula 5.10\n",
    "            g_j = y_j_cap * (np.ones_like(y_j_cap)-y_j_cap)*(y_j - y_j_cap)\n",
    "            # formula 5.15\n",
    "            e_h = b_h *(np.ones_like(b_h)-b_h) * np.dot(g_j, w_h_j.T)\n",
    "\n",
    "            delta_w_h_j = eta * np.dot(b_h.T, g_j)\n",
    "            delta_theta_j = -eta * g_j\n",
    "            delta_v_i_h = eta * (np.dot(e_h.T, x_i)).T\n",
    "            delta_gamma_h = -eta * e_h\n",
    "\n",
    "            w_h_j += delta_w_h_j\n",
    "            theta_j += delta_theta_j\n",
    "            v_i_h += delta_v_i_h\n",
    "            gamma_h += delta_gamma_h\n",
    "\n",
    "            error += 0.5*np.sum((y_j_cap - y_j)*(y_j_cap - y_j))\n",
    "        if i % 10 == 0:\n",
    "            error_list.append(error/len(dataset_array))   \n",
    "    print(error_list)\n",
    "    return w_h_j, v_i_h, theta_j, gamma_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.21627199610317085, 0.1847302176018347, 0.17494199178601993, 0.16656512135352605, 0.16316608869214269, 0.16196825401081952, 0.16138577931454817, 0.16095433576224158, 0.16055459951763884, 0.16016020456760566]\n",
      "0.8120805369127517\n",
      "Wall time: 3.01 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "test_num = dataset.shape[0]//5\n",
    "w_h_j, v_i_h, theta_j, gamma_h = standard_bp(dataset_array[test_num:], niter =100, eta=0.1)\n",
    "acc = cal_acc(dataset_array[:test_num], w_h_j, v_i_h, theta_j, gamma_h)\n",
    "print(acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
